home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 10 / AACD 10.iso / AACD / Online / SpeakFreely / src / lpc10 / invert.c < prev    next >
C/C++ Source or Header  |  2000-05-18  |  2KB  |  76 lines
  1. /*****************************************************************
  2. *
  3. *    INVERT Version 45G
  4. *
  5. *****************************************************************
  6. *
  7. *  Invert a covariance matrix using Choleski decomposition method
  8. *
  9. *  Inputs:
  10. *    ORDER          - Analysis ORDER
  11. *    PHI(ORDER,ORDER) - Covariance matrix
  12. *    PSI(ORDER)       - Column vector to be predicted
  13. *  Outputs:
  14. *    RC(ORDER)          - Pseudo reflection coefficients
  15. *  Internal:
  16. *    V(ORDER,ORDER)   - Temporary matrix
  17. *    X(ORDER)          - Column scaling factors
  18. *
  19. *  NOTE: Temporary matrix V is not needed and may be replaced
  20. *    by PHI if the original PHI values do not need to be preserved.
  21. */
  22.  
  23. #include "config.ch"
  24. #include "lpcdefs.h"
  25. #include <math.h>
  26.  
  27.  
  28. invert(phi, psi, rc )
  29. float phi[MAXORD][MAXORD], psi[], rc[MAXORD][AF];
  30. {
  31. int i, j, k;
  32. static float save, eps=1.0e-10;
  33.  
  34. /*  Decompose PHI into V * D * V' where V is a triangular matrix whose
  35. *   main diagonal elements are all 1, V' is the transpose of V, and
  36. *   D is a vector.  Here D(n) is stored in location V(n,n).        */
  37.  
  38.  
  39. for(j=0;j<ORDER;j++)    {
  40. /**  for(i=j;i<ORDER;i++) {
  41.         v[i][j] = phi[i][j];
  42.     }**/
  43.     for(k=0;k<j;k++)    {
  44.         save = phi[j][k]*phi[k][k];
  45.         for(i=j;i<ORDER;i++)    {
  46.     phi[i][j] -= phi[i][k]*save;
  47.         }
  48.     }
  49.     
  50. /*  Compute intermediate results, which are similar to RC's     */
  51.  
  52.     if((float) fabs((double)phi[j][j]) < eps ) break;
  53.   rc[j][AF-1] = psi[j];
  54.     for(k=0;k<j;k++)
  55.     {
  56.         rc[j][AF-1] -= rc[k][AF-1]*phi[j][k];
  57.     }
  58.     phi[j][j] = 1./phi[j][j];
  59.     rc[j][AF-1] = rc[j][AF-1]*phi[j][j];
  60.     rc[j][AF-1] = mmax(mmin(rc[j][AF-1],.999),-.999);
  61.     
  62.     
  63. }
  64. if((float) fabs((double)phi[j][j]) < eps )    {
  65.  
  66. /*  Zero out higher order RC's if algorithm terminated early    */
  67.  
  68.     for(i=j;i<ORDER;i++)
  69.        rc[i][AF-1] = 0.;
  70.  
  71. }
  72.  
  73.  
  74.  
  75. }
  76.